using System;
using System.Collections.Generic;
using System.IO;
using System.Text;
using System.Xml;

using Atomic.Interfaces.Yphon;
using Atomic.Libraries;
using Atomic.Libraries.Mathematics;
using Atomic.Structures;

namespace Atomic.Vasp.LinearResponse
{
	[Serializable]
	public class VaspLinearResponseResult : VaspResult, IVaspLinearResponseResult, IYphonVaspLinearResponseResult
	{
		[NonSerialized] private bool parsedLinearResponse;
		[NonSerialized] private SiteMatrix forceConstants;

		public VaspLinearResponseResult(DirectoryInfo directory)
			: base(directory)
		{
		}

		public void ParseLinearResponse()
		{
			if (!parsedLinearResponse)
			{
				parsedLinearResponse = true;

				// Must do the basic parsing first.
				Parse();

				XmlDocument xml = new XmlDocument();
				xml.LoadXml(Encoding.Default.GetString(Vasprun.GetData()));

				ReadVasprun(xml);
			}
		}

		private void ReadVasprun(XmlDocument xml)
		{
			Parse(xml["modeling"]);
		}

		private void Parse(XmlNode modelingNode)
		{
			int n = 0;
			List<double> masses = new List<double>();
			foreach (XmlNode atomNode in modelingNode["atominfo"].FindNode("array", "atomtypes")["set"])
			{
				if (atomNode.Name == "rc")
				{
					for (int i = 0; i < int.Parse(atomNode.ChildNodes[0].InnerText); i++)
					{
						masses.Add(atomNode.ChildNodes[2].ToDouble());
						n++;
					}
				}
			}

			if (n != FinalStructure.Sites.Count)
			{
				throw new Exception("Inconsistent number of atoms.");
			}

			// Find last instance of the Hessian matrix.
			XmlNode hessianNode = null;
			foreach (XmlNode node in modelingNode)
			{
				if (node.Name == "calculation" && node["dynmat"] != null && node["dynmat"].FindNode("varray", "hessian") != null)
				{
					hessianNode = node["dynmat"].FindNode("varray", "hessian");
				}
			}

			Matrix a = Matrix.Zero(3 * n, 3 * n);
			for (int i = 0; i < 3 * n; i++)
			{
				a = a.SetMatrix(i, 0, Matrix.Transpose(hessianNode.ChildNodes[i].ToVector()));
			}

			SpaceMatrix[,] forceConstants = new SpaceMatrix[n, n];
			for (int i = 0; i < n; i++)
			{
				for (int j = 0; j < n; j++)
				{
					// Extract matrix of force constants for this pair of atoms; 3x3 submatrix corresponding to the xyz coordinates.
					forceConstants[i, j] = new SpaceMatrix(a.GetMatrix(3 * i, 3 * j, 3, 3).ToArray());

					// Inverse normalization by atomic weights.
					forceConstants[i, j] *= -Math.Sqrt(masses[i] * masses[j]);
				}
			}

			this.forceConstants = new SiteMatrix(FinalStructure.Sites, forceConstants);
		}

		/*/// <summary>
		/// The total elastic tensor. Follows the usual engineering convention.
		/// http://en.wikipedia.org/w/index.php?title=Linear_elasticity&oldid=573528883#Anisotropic_homogeneous_media
		/// http://physics.unl.edu/tsymbal/teaching/SSP-927/Section%2004_Elastic_Properties.pdf
		/// </summary>
		public Matrix ElasticTensor
		{
			get
			{
				throw new NotImplementedException();
			}
		}*/

		/// <summary>
		/// The Hessian of the potential energy. Not normalized by atomic masses.
		/// </summary>
		public SiteMatrix ForceConstants
		{
			get
			{
				ParseLinearResponse();
				return forceConstants;
			}
		}

		IBinaryData IYphonVaspLinearResponseResult.Contcar
		{
			get
			{
				return Contcar;
			}
		}

		IBinaryData IYphonVaspLinearResponseResult.Outcar
		{
			get
			{
				return Outcar;
			}
		}

		/*public LinkedMatrix<Site> PotentialEnergyHessian
		{
			get;
			private set;
		}

		public LinkedVector<Site> Eigenvalues
		{
			get;
			private set;
		}

		public LinkedMatrix<Site> Eigenvectors
		{
			get;
			private set;
		}*/
	}
}
